Read in Data for Myeloarchitecture Characterization

Final study sample

participants <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/sample_info/7T_MP2RAGE_finalsample_demographics.csv") #215 scans from 140 subjects. csv generated by /sample_construction/finalsample_7Tmyelin.Rmd

Glasser regions list

glasser.regions <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist.csv")
glasser.snr.exclude <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/SNR/glasser_SNR_exclusion.csv")
glasser.frontal <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist_frontallobe.csv")
glasser.plotting <- glasser.regions
glasser.plotting$cortex <- "cortex"
glasser.plotting$cortex[glasser.plotting$orig_parcelname %in% glasser.snr.exclude$orig_parcelname] <- "exclude"
glasser.plotting <- glasser.plotting %>% select(label, cortex)
glasser.plotting <- rbind(glasser.plotting, data.frame(label = "rh_???", cortex = "exclude"))
glasser.plotting <- rbind(glasser.plotting, data.frame(label = "lh_???", cortex = "exclude"))
glasser.plotting <- glasser.plotting %>% filter(label != "lh_L_10pp")

Depth-dependent R1 in HCP-MMP regions

#R1 by region matrices at 7 cortical depths
myelin.glasser.7T <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/depthR1_glasseratlas_finalsample.RDS") #generated by /surface_metrics/surface_measures/extract_depthdependent_R1.R
ordered_depths <- c("depth_7", "depth_6", "depth_5", "depth_4", "depth_3", "depth_2", "depth_1")
depth_colorbar <- c("#004A38FF", "#006458FF", "#0093B0FF", "#8AABD5FF", "#C0BFE3FF", "#E0D5ECFF", "#F2E7F4FF")

Regional average R1 in HCP-MMP regions

#Calculate across-depths average regional R1
regionalaverage.myelin <-  do.call(rbind, myelin.glasser.7T) #R1 at all depths
regionalaverage.myelin <- regionalaverage.myelin %>% select(contains("ROI"))
regionalaverage.myelin <- colMeans(regionalaverage.myelin) %>% as.data.frame()  %>% set_names("mean.R1")
regionalaverage.myelin <- regionalaverage.myelin %>% mutate(orig_parcelname = row.names(regionalaverage.myelin))

Myelin Maps

#Myelin Water Fraction Imaging
MWF.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/MyelinWaterFraction_Liu2019/source-Liu2019_desc-MWFmap_space-fsaverage_den-164k_atlas-glasser360.pscalar.nii")
MWF <- data.frame(MWF = MWF.cifti$data, orig_parcelname = names(MWF.cifti$Parcel))

#Myelin Basic Protein Gene Expression
MBP.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/AHBA_magicc/expression_maps/source-magicc_desc-MBPexpression_space-fsLR_den-32k.pscalar.nii")
MBP <- data.frame(MBP = MBP.cifti$data, orig_parcelname = names(MBP.cifti$Parcel))

#Sensorimotor-Association Axis
SAaxis.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii")
SAaxis <- data.frame(SA.axis = rank(SAaxis.cifti$data), orig_parcelname = names(SAaxis.cifti$Parcel))
df.list <- list(regionalaverage.myelin, MWF, MBP, SAaxis) #dfs to merge
myelin.maps <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list) 
myelin.maps <- merge(myelin.maps, glasser.regions, by = c("orig_parcelname"), sort = F)
myelin.maps <- myelin.maps[!(myelin.maps$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels

write.csv(myelin.maps, "/Volumes/Hera/Projects/corticalmyelin_development/code/corticalmyelin_maturation/results/R1_corticalmyelin_anatomy/myelin_maps/myelin_maps.csv", quote = F, row.names = F)

R1 Indexes Cortical Myelin Content Across Cortical Regions

myelin.colorbar <- c("#edf6ff", "#c2d4ed", "#809dc4", "#345c91", "#0e4669")

Whole-brain Map of R1

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = mean.R1), , colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.54, 0.58), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/R1_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

R1 Correlation with Myelin Water Fraction Imaging

MWF Map

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MWF), colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.02, 0.05), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MWF_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

Correlation

cor.test(myelin.maps$mean.R1, myelin.maps$MWF, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  myelin.maps$mean.R1 and myelin.maps$MWF
## S = 2261576, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.6422757

Correlation Plot

myelin.maps %>% filter(MWF < 0.08) %>% #remove 2 MWF major outliers in PHA1
ggplot(., aes(x = MWF, y = mean.R1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nMyelin Water Fraction", y="R1\n") +
  scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MWF_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

R1 Correlation with Myelin Basic Protein Gene Expression

MBP Map

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MBP), colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, limits = c(-0.65, 1), oob = squish, na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MBP_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

Correlation

cor.test(myelin.maps$mean.R1, myelin.maps$MBP, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  myelin.maps$mean.R1 and myelin.maps$MBP
## S = 379704, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5279897

Correlation Plot

myelin.maps %>% 
ggplot(., aes(x = MBP, y = mean.R1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nMyelin Basic Protein Expression (z)", y="R1\n") +
  scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MBP_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

R1 Correlation with the Sensorimotor-Association Axis

S-A axis

myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = SA.axis), colour=I("gray50"), size=I(.06)) + 
scale_fill_gradientn(colors = myelin.colorbar, trans = 'reverse', na.value = "gray75") + theme(legend.position = "none")

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/SAaxis_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)

Correlation

cor.test(myelin.maps$mean.R1, myelin.maps$SA.axis, method = c("spearman"))
## 
##  Spearman's rank correlation rho
## 
## data:  myelin.maps$mean.R1 and myelin.maps$SA.axis
## S = 10280562, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.6261257

Correlation Plot

myelin.maps %>%
ggplot(., aes(x = SA.axis, y = mean.R1)) +
  geom_point(color = c("#8AABD5FF"), size = 0.2) +
  geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
  theme_classic() +
  labs(x="\nSensorimotor-Association Axis", y="R1\n") +
  scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
  theme(
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/SAaxis_R1_correlationplot.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

R1 Indexes Cortical Myelin Content Across Cortical Depths

Studying Cortical Depths Dominated by Signal from Cortical Gray Matter

#Compute average frontal cortex volume fraction at each depth for final study participants
cortexpve <- read_parquet("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/7T_BrainMechanisms_cortexpve.parquet") #parquet with cortex fraction measures
cortexpve <- cortexpve %>% filter(atlas == "glasser") #just for glasser atlas 
cortexpve <- cortexpve[!(cortexpve$StructName %in% glasser.snr.exclude$orig_parcelname),] #just for good regions
cortexpve <- cortexpve[cortexpve$StructName %in% glasser.frontal$orig_parcelname,] #just for frontal lobe

extract_depth_cortexpve <- function(mymeasure){
  measure.df <- cortexpve %>% select(subject_id, session_id, StructName, all_of(mymeasure))
  names(measure.df) <- c("subject_id", "session_id", "orig_parcelname", "cortex.fraction")
  depth.cortexpve <- merge(measure.df, participants, by = c("subject_id", "session_id")) %>% select(cortex.fraction) %>% colMeans(.$cortex.fraction) #get cortex.fraction for final study participants only, and then average for this depth
  return(depth.cortexpve)}

pve.measures <- list("Mean_cortexpve.1.00%", "Mean_cortexpve.0.9%", "Mean_cortexpve.0.8%","Mean_cortexpve.0.7%","Mean_cortexpve.0.6%","Mean_cortexpve.0.5%","Mean_cortexpve.0.4%","Mean_cortexpve.0.3%","Mean_cortexpve.0.2%","Mean_cortexpve.0.1%", "Mean_cortexpve.0.0%") #measures to extract data for 
cortexpve.frontallobe <- lapply(pve.measures, function(x) { 
  extract_depth_cortexpve(x)}) 
cortexpve.frontallobe <- do.call(rbind, cortexpve.frontallobe) %>% as.data.frame()
cortexpve.frontallobe$depth <- c("a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k")
cortexpve.frontallobe$depth <- factor(cortexpve.frontallobe$depth, ordered = T, levels = c("k", "j", "i", "h", "g", "f", "e", "d", "c", "b", "a"))
kable(cortexpve.frontallobe)
cortex.fraction depth
0.7103157 a
0.8388772 b
0.9281964 c
0.9749537 d
0.9922900 e
0.9960812 f
0.9925421 g
0.9710414 h
0.9071242 i
0.7684002 j
0.5157232 k
cortexpve.frontallobe %>% 
  ggplot(aes(x = cortex.fraction, y = depth)) +
  geom_point(color = "#809dc4", shape = 15, size = 2) +
  theme_classic() +
  xlab("\nCortex Fraction") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text.x = element_text(size = 6, family = "Arial", color = c("black")),
  axis.text.y = element_blank(),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks.x = element_line(linewidth = .2),
  axis.ticks.y = element_blank()) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1_supplementary/Depthplot_cortexfraction_frontal.pdf", device = "pdf", dpi = 500, width = 2.5, height = 2.5)

Regional R1 Distributions by Depth

#Compute mean R1 in each region at each cortical depth
regional_means <- function(input.df){
  meanmap <- input.df %>% select(contains("ROI")) %>% colMeans() %>% as.data.frame() %>% set_names("mean.R1")
  meanmap$orig_parcelname <- row.names(meanmap)
  meanmap <- merge(meanmap, glasser.regions, by = c("orig_parcelname"), sort = F)
  return(meanmap)
}

regional.myelin.alldepths <- lapply(myelin.glasser.7T, function(depth) {
  regional_means(depth)
})
regional.myelin.alldepths <-  do.call(rbind, regional.myelin.alldepths) #regional mean R1 at all depths
regional.myelin.alldepths$depth <- substr(row.names(regional.myelin.alldepths), 1, 7) #assign depth variable
regional.myelin.alldepths$depth <- factor(regional.myelin.alldepths$depth, ordered = T, levels = ordered_depths)

regional.myelin.alldepths <-  regional.myelin.alldepths[!(regional.myelin.alldepths$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels
regional.myelin.alldepths <- regional.myelin.alldepths[regional.myelin.alldepths$orig_parcelname %in% glasser.frontal$orig_parcelname,] #only studying frontal


regional.myelin.alldepths %>% 
  ggplot(aes(x = mean.R1, y = depth)) +
  geom_violin(aes(fill = depth), color = "white") +
  stat_summary(fun = median, geom = "point", shape = 23, size = 2.25, color = "white") +
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  scale_x_continuous(limits = c(0.49, 0.65), expand = c(0,0)) +
  scale_y_discrete(expand = c(0, 0)) +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Depthplot_regionalR1_frontal.pdf", device = "pdf", dpi = 500, width = 2.6, height = 4.9)
for(this.depth in c(unique(regional.myelin.alldepths$depth))){
  
  depth.data <- regional.myelin.alldepths %>% filter(depth == this.depth) %>% filter(label != "lh_L_10pp")
  
  depth.plot <- ggplot() + 
  geom_brain(data = depth.data, atlas = glasser,  mapping = aes(fill = mean.R1), colour=I("gray50"), size=I(.06)) + theme_classic() + theme(legend.position = "none") + 
   paletteer::scale_fill_paletteer_c("grDevices::PuBuGn", direction = -1, limits = c(0.5, 0.62), oob = squish, na.value = "white") + theme(legend.position = "none") +
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray50", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75")) 

  print(depth.plot)
  ggsave(filename = sprintf("/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/%s_corticalmap_frontal.pdf", this.depth), device = "pdf", dpi = 500, width = 4.95, height = 2.25)
}

R1 Myelin Profile in Individual Regions

Group Average Profiles

examplar.areas.alldepth <- rbind(area4.myelin.alldepths, area46.myelin.alldepths, areaa24.myelin.alldepths)

ggplot() +
  geom_smooth(data = examplar.areas.alldepth, aes(x = mean.R1, y = depth, group = orig_parcelname, color= orig_parcelname), se = F, linewidth = .8) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  theme_classic() +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Area_profiles_group.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

Individual Profiles

#R1 values in all regions at all depths for all participants
myelin.glasser.7T.alldepths <-  do.call(rbind, myelin.glasser.7T)
myelin.glasser.7T.alldepths$depth <- substr(row.names(myelin.glasser.7T.alldepths), 1, 7) #assign depth variable
myelin.glasser.7T.alldepths$depth <- factor(myelin.glasser.7T.alldepths$depth, ordered = T, levels = ordered_depths)

#Wide to long format
cols_to_pivot <- names(myelin.glasser.7T.alldepths)[grepl("ROI", names(myelin.glasser.7T.alldepths))]
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% pivot_longer(cols = all_of(cols_to_pivot), names_to = "orig_parcelname", values_to = "R1")

#Select exemplar regions of interest
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI" | orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI" | orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI")

#Format grouping variables
myelin.glasser.7T.alldepths$subject_id <- as.factor(myelin.glasser.7T.alldepths$subject_id)
myelin.glasser.7T.alldepths$orig_parcelname <- as.factor(myelin.glasser.7T.alldepths$orig_parcelname)

#Plot for the oldest 3 participants 
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(age > 32)
ggplot() +
  geom_smooth(data = myelin.glasser.7T.alldepths, aes(x = R1, y = depth, color = orig_parcelname, group = interaction(orig_parcelname, subject_id)), se = F, linewidth = .4) + 
  scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
  theme_classic() +
  xlab("\nR1") +
  ylab("Cortical Depth\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Area_profiles_individual.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)

Region-wise R1 Correlations Between Depths

regional.myelin.alldepths.wide <- regional.myelin.alldepths %>% pivot_wider(id_cols = orig_parcelname, names_from = depth, values_from = mean.R1) %>% select(-orig_parcelname)
regional.myelin.alldepths.wide <- regional.myelin.alldepths.wide %>% select(depth_7, depth_6, depth_5, depth_4, depth_3, depth_2, depth_1)
depth.R1.cors <- cor(regional.myelin.alldepths.wide) # compute correlation matrix

ggcorrplot(corr = depth.R1.cors, method = "square", type = "upper", show.diag = T, lab = F, lab_size = ) +
scale_fill_gradientn(colors = c("#edf6ff", "#c2d4ed", "#97adcc", "#345c91"))

ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Depth_correlationmatrix_frontal.pdf", device = "pdf", dpi = 500, width = 4.4, height = 2.2)